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for  this  investigation  and  as  the  source  of  a  case  study  we  describe  in  the 
paper.  The  principal  focus  of  the  paper  is  on  the  development  of  an  approach 
that  overcomes  the  combinatorial  explosion  of  truly  optimal  estimation 
algorithms.  We  accomplish  this  by  constructing  a  systematic  design 
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that  we  address-  concern  the  way  in  which  these  estimators  interact  and  the 
method  each  estimator  uses  to  account  in  its  own  model  for  the  influence  of 
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I.  Introduction 


In  Doerschuk  (1986)  we  developed  a  methodology  for  modeling 
electrocardiograms  (ECG’s)  that  could  be  used  as  the  basis  for  EOG  signal 
processing/analysis  algorithms.  The  models  we  have  developed  have  several 
important  characteristics: 


(1)  The  models  are  hierarchical  in  nature.  Specifically,  at  the 
upper  level  we  model  the  event  structure  of  the  heart,  capturing 
the  electrical  state  of  various  anatomical  portions  of  the  heart. 
At  the  lower  level  we  model  the  actual  observation,  the  EOG 
signal.  The  model  reflects  the  fact  that  particular  cardiac 
events  directly  result  in  the  various  waveforms  seen  in  the  EOG 
by  having  particular  changes  of  state  in  the  upper  model  initiate 
waveforms  in  the  lower  level  model. 

(2)  The  upper  level  model  consists  of  interacting  finite-state 
processes,  each  of  which  models  a  specific  anatomical  portion  of 
the  heart.  In  this  way  we  attempt  to  capture  the  distributed  but 
coordinated  way  in  which  the  heart  operates.  In  particular,  this 
model  structure  allows  us  to  highlight  the  aspects  of  timing  and 
control  that  are  critical  to  cardiac  behavior.  Specifically,  in 
our  models  the  state  transition  probabilities  of  each  subprocess 
are  affected  by  the  states  of  other  subprocesses,  allowing  us  to 
model,  for  example,  the  attempt  of  one  portion  of  the  heart  to 
precipitate  sin  event  (e.g.  a  contraction)  in  smother  portion  of 
the  heart. 

For  a  detailed  description  of  this  modeling  methodology  see  [Doerschuk  1985, 


1986].  An  example  is  also  described  in  the  next  section  as  we  develop  our 
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approach  to  estimation. 

The  motivation  for  developing  models  of  this  type  was  to  overcome 
limitations  of  existing  signal  processing  models  (e.g.  those  used  by  Gustafson 
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1978a.  b.  1981;  Ciocloda  1983;  Gersch  1970,  1975;  Tsui  1975;  Grove  1978; 
Haywood  1977;  Richardson  1976)  by  capturing  physiology  in  a  more  fundamental 
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way  while,  however,  stopping  far  short  of  the  level  of  detail  found  in 
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physiologically  accurate  models.  We  have  attempted  to  develop  models  that 
capture  all  features  needed  to  characterize  differences  in  cardiac  rhythms, 
i.e.  differences  in  the  sequential  nature  of  various  cardiac  events. 
Consequently,  the  focus  of  our  Interest  in  this  paper  is  on  the  estimation  of 
the  states  of  the  upper,  level  of  the  model.  As  discussed  in  [Doerschuk  1986], 
previously  developed  rhythm  analysis  methods  have  also  been  based  on 
event-level  models  of  cardiac  behavior.  These  models  differ  from  ours  in 
several  Important  respects.  Specifically,  in  previous  methods  typically  only 
the  event  level  description  is  modeled,  and  it  is  assumed  that  event-level 
Inputs  are  available  from  a  wave  detection  preprocessor.  For  methods  in  which 
attention  is  paid  only  to  ventricular  events  —  i.e.  to  so-called  R-waves  — 
extremely  useful  models  of  this  type  can  be  developed  since  (a)  R-waves  are 
comparatively  large  and  thus  can  be  detected  with  great  reliability;  and  (b) 
one  can  achieve  great  efficiencies  in  describing  event-level  models  by  using 
as  a  sequential  index  the  successive  R-wave  occurrences  rather  than  the  real 
EOG  sample-to-sample  time.  However,  when  one  begins  to  include  more  cardiac 
detail,  such  as  the  P-waves  arising  from  atrial  activity,  difficulties  arise 
with  such  high-level  models.  In  particular,  because  of  their  much  lower 
energy  levels  P-waves  are  much  more  difficult  to  detect.  Consequently,  in 
some  previous  methods  one  finds  ad  hoc  attempts  at  feedback  from  the  rhythm 
tracker  to  the  wave  detection  process,  taking  advantage  of  the  fact  that  in  a 
normal  heartbeat  the  P-wave  precedes  the  R-wave  in  a  predictable  fashion. 
Furthermore,  it  is  very  difficult  in  such  a  framework  to  incorporate  the 
possibility  of  missed  detections  or  false  alarms,  which  are  very  real 
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possibilities,  especially  for  P-waves  and  especially  in  aberrant  cardiac 
rhythms.  For  our  two-level  models  in  which  the  measurements  are  the  raw  data, 
wave  detection  feedback  and  the  possibilities  of  false  positive  and  missed 
detections  are  built  in  in  a  fundamental  way.  Also,  while  in  a  perfectly 
normal  rhythm  P-waves  and  R-waves  are  always  paired,  cardiac  arrhythmias  are 
characterized  by  asynchronous  behavior,  in  which  for  example  several  P-waves 
may  occur  without  an  intervening  R-wave  or  vice  versa.  The  emphasis  in  our 
models  on  timing  captures  such  asynchronous  behavior  in  a  natural  way  rather 
than  in  the  less  than  completely  satisfactory  manner  found  in  previous 
methods. 

The  various  potential  advantages  we  have  attributed  to  using  the  models 
in  [Doerschuk  1986]  for  EOG  rhythm  processing  do  not  come,  of  course,  without 
a  significant  price.  In  particular,  while  we  do  contend  that  truly  optimal 
estimation  based  on  these  models  would  achieve  these  advantages,  the 
computational  load  associated  with  optimal  processing  is  prohibitively  large. 
Thus  the  major  issue  is  the  development  of  feasible,  suboptlmal  estimation 
algorithms.  In  this  paper  we  investigate  the  development  of  such  algorithms 
that  take  advantage  of  two  important  features  of  this  class  of  estimation 
problems.  First,  the  estimation  of  event  sequences  in  the  upper  level  model 
is  essentially  a  decoding  problem  (i.e.  the  EOG  is  an  encoding  of  the  discrete 
cardiac  events  we  wish  to  estimate).  Consequently  we  make  repeated  use  of  an 
efficient  technique  for  optimal  estimation  of  finite-state  processes  first 
developed  for  coding  applications,  namely  the  Viterbi  algorithm  [Forney  1973]. 
Second,  since  our  models  are  distributed,  we  can  consider  the  design  of 
distributed  estimators,  consisting  of  interacting  algorithms  each  focused  on 
the  job  of  estimating  the  state  of  a  particular  subprocess.  Such  estimation 
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structures  offer  the  attractive  possibility  of  implementation  in  a  distributed 
processor,  thereby  allowing  significant  improvements  in  throughput  rates. 

The  design  of  such  estimators  also  raises  a  number  of  important 
questions.  In  particular,  since  the  several  subprocesses  of  our  upper  level 
model  interact  strongly,  it  is  not  possible  to  estimate  the  state  of  a 
'  .^process  without  accounting  for  the  influence  on  it  of  other  subprocesses. 
Consequently  it  is  necessary  to  include  an  (hopefully  highly)  aggregated  model 
of  other  subprocesses  that  captures  the  dynamics  of  the  interactions  these 
subprocesses  have  with  the  particular  subprocess  being  estimated.  Also,  it  is 
necessary  for  the  estimators  of  interacting  subprocesses  to  interact 
themselves  (e.g.  estimators  of  atrial  and  ventricular  activity  most  certainly 
have  information  worth  sharing!).  This  raises  an  interesting  problem. 
Specifically,  since  each  estimator  in  essence  receives  additional  measurements 


in  the  form  of  information  passed  from  other  estimators,  it  is  necessary  to 
model  these  measurements  —  l.e.  each  estimator  needs  an  aggregated  model  of 
the  dynamics  and  the  uncertainties  in  the  other  estimators.  In  addition, 
since  each  estimator  is  using  the  same  raw  data  (the  EOG)  but  is  interested  in 
only  some  of  the  events  in  the  data,  it  may  be  necessary  to  provide 
Information  to  each  estimator  concerning  estimated  times  of  occurrence  of 
other  events  in  the  EOG  data  (e.g.  an  atrial  estimator  may  need  estimates  of 
R-wave  locations  from  the  ventricular  estimator  in  order  to  assist  it  in 
locating  the  much  smaller  P-waves).  Also,  as  one  might  expect,  there  may  very 
well  be  a  need  for  some  iteration  in  this  process  so  that  a  high  level  of 
performance  and  consistency  among  the  estimators  is  achieved. 

While  electrocardiogram  analysis  has  provided  the  motivation  and  examples 


for  our  work,  there  are  a  variety  of  other  applications  in  which  similar 
estimation  problems  arise.  In  particular,  consider  the  distributed  monitoring 
of  Interconnected  power  systems.  Such  systems  are  made  up  of  strongly 
Interacting  components  subject  to  events  (such  as  generator  trips  and  line 
faults)  that  can  precipitate  events  in  other  parts  of  the  system.  An 
extremely  important  problem  is  the  design  of  distributed  monitoring  systems, 
and  a  critical  aspect  of  this  problem  is  determining  how  to  structure  the 
interaction  among  local  monitoring  systems  in  order  to  produce  a  consistent 
and  accurate  overall  estimate  of  system  status.  Similar  issues  also  arise  in 
military  contexts  in  distributed  battle  management  and  assessment.  In  all  of 
these  problems  the  key  question  is  how  do  we  design  distributed  event 
estimation  algorithms,  and  in  this  paper  we  address  this  problem  and  in 
particular  the  several  critical  issues  raised  in  the  previous  paragraph.  We 
begin  in  the  next  section  with  a  case  study,  in  which  we  describe  the  details 
of  designing  an  estimator  for  an  example  corresponding  to  a  particular  cardiac 
arrhythmia.  This  case  study  allows  us  to  introduce  the  major  questions  that 
arise  in  designing  distributed  event  estimation  algorithms.  In  Section  3  we 
then  step  back  from  this  example  and  extract  from  it  a  general,  systematic 
design  approach.  We  then  discuss  a  series  of  simplifications  that  are 
possible  for  the  class  of  models  arising  in  EOG  analysis  and  that  result  in 
the  algorithm  described  in  Section  2. 

While  the  motivation  for  considering  distributed  estimators  for  EOG 
rhythm  analysis  is  the  issue  of  computational  complexity  (as  opposed  to 
geographic  separation  of  sensors  and  controls  as  in  the  case  of  electric  power 
or  military  systems),  we  have  not  taken  our  design  methodology  to  the  point  of 


developing  a  complete  and  computationally  feasible  EOG  analysis  system,  and 
indeed  as  we  discuss  in  Section  4  a  number  of  issues  remain  before  such  a 
system  is  developed.  Rather  what  we  have  attempted  to  do  is  to  use  the  EOG 
problem  as  a  context  in  which  we  can  explore  the  major  problems  in  the  design 
of  distributed  event  estimation  systems,  present  the  elements  needed  in  their 
solution  and  a  systematic  procedure  for  their  application,  and  demonstrate  the 
potential  of  this  approach  for  solving  a  class  of  extremely  complex  estimation 
problems. 

2.  An  Estimation  Example 

In  this  section  we  Introduce  our  approach  to  estimation  algorithm  design 
by  considering  a  particular  example.  The  process  whose  state  is  to  be 
estimated  and  which  is  illustrated  in  Figure  1  models  normal  cardiac  rhythm 
with  occasional  reentrant-mechanism  premature  ventricular  contractions  (PVCs; 
these  result  from  a  normal  excitation  of  the  ventricles  in  effect  circling 
back  on  itself  and  causing  additional  ventricular  contractions).  Let  us 
briefly  indicate  several  important  features  of  the  model: 

(1)  The  model  consists  of  two  subprocesses  one  (the  SA-atrial 

submodel,  denoted  CO,  with  state  x^)  representing  the  behavior  of 

the  upper  chambers  of  the  heart  and  the  other  (the  AV-ventricular 
submodel,  denoted  Cl,  with  state  capturing  the  behavior  of 

the  atrial-ventricular  connection  and  the  lower  chambers  of  the 
heart.  In  the  figure  we  have  depicted  the  state  transition 
probabilities  for  each  submodel  and  the  interaction  between 
submodels  which  is  captured  by  the  dependence  of  these 
probabilities  on  the  state  of  the  other  subprocess.  We  have  also 
indicated  which  transitions  in  each  subprocess  give  rise  to 
waveforms  or  signatures  that  appear  in  the  observed  EOG.  The 
signatures  modeled  are  the  P-wave  (labeled  Pq,  Pj ,  P^.  P^  in  the 
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Figure  1:  A  model  of  normal  cardiac  rhythm  with  occasional  reentrant- 
mechanism  PVC's:  (a)  the  two  subprocesses ; 


figure),  the  R-wave ,  the  T-wave  (corresponding  to  ventricular 
repolarization  following  a  normal  contraction),  and  a  V-wave 
corresponding  to  an  aberrant  reentrant  PVC.  The  actual  signature 
that  appears  in  the  EOS  consists  of  a  mean  value  plus  additive 
white  noise  of  specified  variance.  Each  occurrence  of  these 
signatures  Includes  a  different  realization  of  the  noise.  The 
means  and  variances  for  the  various  signatures  are  also 
illustrated  in  Figure  1. 

(2)  The  interactions  between  the  submodels  are  infrequent  but  are 
extremely  strong.  In  particular,  the  diagram  shown  for  the 
SA-atrial  submodel  represents  normal  activity  which  occurs  unless 
Xj  =  13  in  the  AV-ventricular  submodel,  corresponding  to  the 

initiation  of  a  PVC.  When  such  an  event  occurs,  it  is  possible 
for  the  electrical  signal  to  propagate  back  to  the  upper  chambers 
of  the  heart  and  in  essence  reset  the  timing  of  the  heart's  own 
pacemaker.  This  is  captured  by  modifying  the  transition 
probabilities  of  Xq  so  that  with  probability  1/2  xQ  is  reset  to 

state  25  when  x^=  13,  and  with  probability  1/2  proceeds  in  a 

normal  fashion.  In  the  x^  submodel  the  only  transition 

probability  affected  by  the  value  of  x^  is  p^j.  In  particular, 

Xj  =  0  represents  the  resting  state  of  the  ventricles,  which  is  a 

trapping  state  (p^  =  0)  until  the  ventricles  are  excited  by  an 
atrial  contraction.  This  event  is  modeled  by  x^  =  0  which  both 
initiates  the  P-wave  and  excites  the  AV-ventricular  submodel  by 
causing  p^  to  change  to  a  value  of  1  for  one  time  step. 

(3)  In  our  model,  the  EOG  measurements  are  available  at  a  rate  four 
times  the  clock  rate  of  the  Xq,  x^  processes.  In  order  to  make 

it  possible  for  signatures  to  start  at  any  observation  sample, 
each  signature  appears  four  times  with  0,  1,  2,  or  3  leading 
zeros  in  the  mean  and  covariance  sequences.  (The  subscripts  on 
the  wave-types  indicate  the  number  of  leading  zeros).  This 
explains  states  22-24  in  submodel  CO  and  states  22-30  in  submodel 
Cl. 

(4)  The  initiation  of  reentrant  PVC’s  is  modeled  by  transitions  out 
of  states  12  and  21  in  submodel  Cl.  Occupancy  of  state  12 
corresponds  to  the  completion  of  a  normal  R,  T-wave  pair,  and 
from  this  state  there  is  a  probability  of  0.9  of  returning  to  the 
resting  state  and  a  probability  of  0.1  of  entering  state  13 
corresponding  to  the  initiation  of  a  reentrant  PVC.  Note  that 
there  is  a  much  higher  probability  (0.4)  of  initiating 
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subsequent,  consecutive  reentrant  PVC’s  (the  21-to-13  transition) 
which  results  in  occasional  occurrences  of  bursts  of  aberrant 
PVC’s  as  are  seen  in  episodes  of  ventricular  tachycardia. 

(5)  The  remaining  states  and  transition  probabilities  model  cardiac 
timing  —  propagation  delays,  recovery  time  following 
contraction,  etc.  The  model  does  allow  for  some  uncertainty  in 
this  timing  behavior  and  therefore  some  variability  in  the  heart 
rate  (which  with  a  Markov  chain  cycle  time  of  0.04  seconds  is,  on 
the  average,  75  beats  per  minute)  as  can  be  seen  in  the 
transitions  with  probabilities  less  than  1.  It  is  certainly 
possible  to  add  even  more  variability,  but  for  simplicity  we  have 
not  done  that  here. 

Figure  2  shows  a  plot  of  several  typical  segments  of  a  simulated  ECG 
obtained  using  this  model.  Below  the  EOG  tracing  are  several  sets  of 
annotations.  The  top  row  of  annotations  indicates  the  true  times  and  types  of 
waves  that  are  present  in  the  data  (corresponding  to  the  times  at  which 
transitions  are  made  out  of  state  0  in  submodel  CO  (P-wave)  and  states  4 
(R-wave).  7  (T-wave),  and  13  (V-wave)  of  submodel  Cl).  The  remaining  rows 
represent  various  annotations  constructed  during  the  estimation  process,  with 
the  bottom  row  representing  our  final  set  of  estimates.  Before  turning  to  a 
discussion  of  our  estimation  procedure,  let  us  briefly  comment  on  the  nature 
of  the  simulated  EOG  itself.  In  particular,  the  simulated  data  resembles  an 
EOG,  but  it  is  not  entirely  realistic.  The  chief  difference  is  the  character 
of  the  noise.  First,  no  bandlimited  noise,  "shot  noise"  like  events,  or 
baseline  instabilities  are  present.  Second,  the  noise  that  is  present,  though 
it  resembles  electromyogram  artifact  when  considered  locally,  appears  to  have 
a  time-varying  intensity  (higher  during  the  P-,  R-,  T-,  and  V-waves)  which  is 
not  realistic  for  true  electromyogram  artifact.  (As  discussed  in  [Doerschuk 
1986],  the  reason  for  the  time-varying  intensity  was  to  give  a  crude  model  of 
the  morphology  variability  of  the  different  wave  types).  Each  of  these  issues 
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Figure  2:  Several  segments  of  a  simulated  ECG  obtained  using  the  model  in  Figure  1.  Annotations  below 
the  traces  refer  to  estimates  produced  at  several  points  in  the  estimation  algorithm 
(see  text) . 
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can  be  dealt  with  by  adding  additional  detail  to  our  model.  Again  for 
simplicity  we  have  not  done  so  here.  Also,  much  of  this  detail  may  very  well 
not  be  necessary  for  the  purpose  of  designing  rhythm  processing  algorithms. 

Finally,  in  anticipation  of  the  discussion  to  follow,  we  introduce  some 
notation.  Specifically,  it  is  useful  to  have  a  compact  pictorial  description 
of  interacting  Markov  chains.  Such  a  description  is  illustrated  in  Figure  3. 
Here  the  label  00  denotes  the  SA-atrial  submodel  and  Cl  the  AV-ventriclar 
submodel  shown  in  Figure  1.  The  arrows  between  CO  and  Cl  indicate  that  the 


h01(n> 


state  of  each  subprocess  influences  the  transition  behavior  of  the  other. 

Also,  the  arrows  labeled  P,  R.  T.  and  V  indicate  the  waveforms  initiated  by 

each  subprocess*.  In  addition,  we  use  the  variables  hg^n)  to  denote  the 

sequence  of  interactions  initiated  by  00  and  impinging  on  Cl.  That  is  hQ1(n) 

completely  captures  the  Influence  CO  has  on  the  transition  probabilities  of  Cl 

for  the  transition  Xj(n)  *♦  Xj(n+1).  Referring  to  Figure  1.  we  see  that  we  can 

define  hQ1(n)  so  that  it  takes  on  only  two  values 

fO  if  Xg(n)  =  0 

(1) 

otherwise 

The  only  transition  probability  of  Cl  that  is  influenced  by  CO  is  a  trivial 

function  of  hg^(n): 

(1.  bQ1(n)  =  0 

h0i(«)  =  1 

Similarly  we  can  define  the  interactions  h1Q(n)  from  Cl  impinging  on  CO 

!0  if  x. (n)  *  13 

1  otherwise 

so  that  if  h1Q(n)  =  0  the  transition  probabilities  are  as  indicated  in  the 
figure,  and  if  h^g(n)  =  1  they  are  the  average  of  these  values  and  a 
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probability  1  reset  to  state  25  from  any  other  state.  Note  that  there  are  far 
fewer  values  for  these  interaction  variables  than  for  the  corresponding 
states.  This  fact  is  used  in  an  essential  way  in  constructing  several 
aggregate  models  used  in  our  estimation  methodology. 

As  discussed  in  the  Introduction,  our  approach  to  state  estimation  for 
such  a  process  Involves  the  design  of  a  set  of  interacting  estimators,  each  of 
which  focuses  on  estimation  for  a  particular  subprocess.  Also,  as  we 
indicated,  the  existence  of  the  interactions  among  subprocesses  may  require 
some  iteration.  For  the  present  example  our  estimator  can  be  viewed  as 
consisting  of  three  passes: 

(1)  Derive  a  preliminary  estimate  of  ventricular  activity  (submodel 
Cl). 

(2)  Based  on  the  observed  ECG  and  the  estimates  from  pass  1,  compute 
an  estimate  of  atrial  activity  (submodel  CO). 

(3)  Refine  the  ventricular  estimate  based  on  the  observed  ECG  and  the 
estimates  of  atrial  activity  from  pass  2. 

The  results  from  (2)  and  (3)  form  the  final  estimate.  Note  that  this  approach 
roughly  parallels  the  heuristic  approach  humans  take  in  first  identifying  the 
high  signal-to-noise  ratio  (SNR)  events  (R-  and  V-waves),  then  using  these 
estimates  to  assist  in  locating  the  low  SNR  events  (P-waves).  and  finally 
making  adjustments  to  ensure  accuracy  and  consistency.  Also,  while  we 
describe  these  three  steps  as  separate  passes  through  the  data,  it  is 
straightforward  to  construct  a  pipelined  structure  in  which  the  three  steps 
proceed  at  the  same  time. 

We  now  turn  to  a  detailed  examination  of  each  of  these  three  passes. 


Because  the  first  pass  focuses  on  submodel  Cl.  it  is  natural  to  include  an 
exact  copy  of  this  submodel  in  the  estimator’s  model.  However,  it  is  also 
necessary  to  model  the  interactions  impinging  on  the  Cl  submodel,  i.e.  h^^(n). 
Here  one  is  confronted  with  a  range  of  possibilities,  from  the  exact  model  of 
00  depicted  in  Figure  1  (which  provides  complete  accuracy  in  modeling 
interactions  at  the  expense  of  computational  complexity)  to  no  model  (which 
completely  misses  capturing  interactions  but  leads  to  the  simplest  estimator). 
We  have  chosen  in  this  study  to  use  the  simplest  possible  aggregate  model  for 
submodel  CO  with  which  we  can  still  capture  the  full  range  of  interactions 
with  submodel  Cl.  In  particular  the  aggregate  version  of  submodel  CO  that  we 
use  in  this  step  is  a  two-state  model,  corresponding  to  the  two  possible 
values  of  h^^(n).  In  addition,  we  allow  submodel  Cl  to  reset  the  state  of  our 
two-state  aggregate  model,  again  reflecting  behavior  seen  in  the  full  model. 

In  the  full  discussion  of  our  approach  to  estimation,  this  type  of  aggregate 
model  is  referred  to  as  an  "SO-submode 1 " .  Details  for  this  example  are  given 
in  Figure  4. 

There  are  several  further  points  to  make  about  this  first  pass.  Note 
first  of  all  that  the  subprocess  SO  does  not  initiate  P-waves.  While  it  would 
certainly  be  possible  to  include  the  P-wave,  this  wave  has  a  relatively  small 
amplitude  compared  to  the  ventricular  waves  of  primary  concern  in  this  pass. 
Consequently,  there  is  little  chance  of  confusing  an  unmodeled  P-wave  for  tin 
R-  or  V-wave,  and  for  this  reason  we  chose  to  ignore  its  presence  in  this 
pass.  A  second  question  that  arises  is  the  choice  of  the  value  of  p  in 
submodel  SO.  One  can  imagine  several  methods  for  choosing  p  —  matching  some 
statistic  of  the  exact  submodel  00  or  viewing  p  as  a  design  parameter  to  be 
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Figure  4:  Model  for  the  first  pass  of  the  estimation  algorithm?  (a) 
overall  block  diagram;  (b)  detail  of  the  SO  model— -state  0 
corresponds  to  h  fn)  “  0,  state  1  to  h  ( n)  »  1. 


chosen  to  optimize  estimator  performance.  In  [Doerschuk  1985]  several  general 
statistical  methods  (which  can  be  easily  automated)  are  described  for  choosing 
parameters  to  match  particularly  useful  statistics.  In  Section  3  we  describe 
the  statistical  method  used  to  obtain  the  value  for  p  indicated  in  the  figure. 
Finally,  with  this  parameter  specified  we  have  a  complete  model,  and  the  first 
step  estimator  is  designed  to  produced  a  minimum  probability-of-error  state 
trajectory  estimate  for  this  model  (i.e.  estimates  of  the  states  of  SO  and  Cl 
as  functions  of  time)  based  on  the  observed  ECG.  The  method  used  to  compute 
this  sequence  of  estimates  is  the  Viterbi  algorithm  [Forney  1973]  which 
efficiently  and  recursively  eliminates  candidate  trajectories  that  are 
suboptimal.  Note  that  the  estimate  so-produced  is  an  optimal  smoothed 
trajectory,  i.e.  the  best  state  estimate  at  each  time  is  based  on  information 
before  and  after  that  time.  The  Viterbi  algorithm  minimizes  the  extent  of  the 
noncausality  required  in  computing  this  optimal  estimate. 

The  results  of  this  first  pass  estimator  are  illustrated  in  the  second 
row  of  annotations  in  Figure  2,  where  we  have  indicated  the  estimated  times  of 
occurrence  of  R-.  T-.  and  V-waves.  For  the  most  part  these  estimates  are 
quite  accurate,  thanks  to  the  high  SNR  of  these  waves,  although  there  are 
infrequent  false  alarms  in  the  estimates  caused  by  extra-long  P-P  intervals  in 
which  case  the  estimator  attempted  to  match  a  T-wave  with  an  actual  P-wave. 

The  second  step  in  our  overall  estimation  structure  is  to  estimate  the 
state  trajectory  in  the  SA-atrial  submodel.  Therefore,  it  is  natural  to 
Include  an  exact  copy  of  the  SA-atrial  submodel  in  the  estimator’s  model.  For 
estimating  this  trajectory  the  only  direct  information  from  the  EOG  is  the  low 
SNR  P-wave.  However,  there  is  also  a  great  deal  of  indirect  information 


14 


t 


E 


available  through  the  causal  relationship  between  P-  and  R-  waves  and  V-and 

waves . 

First  consider  interactions  initiated  by  submodel  00.  That  is,  consider 
the  causality  between  P-  and  R-waves  the  latter  of  which  only  occur  when  the 
SA-atrlal  submodel  successfully  excites  the  AV-ventrlcular  submodel.  The 
question  in  this  case  is  how  we  can  exploit  the  auxiliary  information 
concerning  R-wave  occurrences  determined  in  the  first  estimation  pass.  At  the 
very  least  one  could  Imagine  using  the  state  estimates  for  SO  which  are 
estimates  of  interactions  impinging  on  the  AV- ventricular  submodel.  Since  the 
O-state  in  this  submodel  corresponds  to  the  0-state  in  the  original  submodel 
CX)  (and  thus  to  attempts  to  excite  submodel  Cl),  the  estimates  of  times  at 
which  SO  is  in  state  0  would  be  likely  estimates  of  times  at  which  h^^n)  =  0. 
However,  because  of  the  highly  aggregated  nature  of  SO.  some  of  these 
estimates  may  be  somewhat  suspect.  However,  when  such  an  estimate  is  coupled 
together  with  a  closely  following  estimated  occurrence  of  an  R-wave 
(corresponding  to  :he  estimate  of  the  Cl  subprocess  occupying  state  4).  the  SO 
estimate  is  much  more  likely  to  correspond  to  a  true  occurrence  of  an  attempt 
at  ventricular  excitation.  Consequently  the  information  we  provide  to  pass  2 
from  pass  1.  which  we  will  refer  to  as  estimated  augmented  interactions, 
consists  of  the  sequence  of  estimates  of  the  states  of  SO  and  Cl  produced  in 
pass  1 . 

We  now  describe  the  issues  involved  in  using  these  estimated  augmented 
interactions.  First  of  all.  these  estimates  contain  errors,  so  some 
uncertainty  in  them  must  be  modeled.  Note,  however,  that  the  errors  of 
importance  here  are  not  only  memoryless  errors  (which  could  be  modeled  by 


static  mlsclassiflcatlon  probabilities)  but  also  errors  in  timing  (e.g.  the 
estimated  time  of  occurrence  of  an  R-wave  may  be  in  error  by  one  or  two 
samples).  Consequently,  we  need  a  dynamic  model  for  the  way  in  which 
estimated  augmented  interactions  provide  information  about  CO.  The  way  in 
which  this  is  accomplished,  as  illustrated  in  Figure  5,  is  by  modeling  the 
estimated  interactions,  denoted  by  z^.  as  the  observed  outputs  of  an 
additional  submodel  of  a  class  we  refer  to  as  SI  submodels.  This  additional 
submodel  receives  interactions  from  submodel  CO.  whose  state  we  wish  to 
estimate.  In  order  to  model  the  fact  that  the  estimates  in  z^(n)  may  contain 
time  shifts  relative  to  the  actual  values  of  the  interactions  h^^n),  we  take 
as  the  state  of  the  SI  submodel  a  vector  of  the  most  recent  interaction 
values.  To  minimize  the  size  of  the  SI  state  space  one  clearly  wishes  to 
minimize  the  dimension  of  this  vector.  For  this  study  we  found  a  dimension  of 
2  to  be  adequate,  so  that  the  state  of  SI  at  time  n  is  (hQ1(n-l).  hQ1(n-2)). 

By  examining  the  CO  submodel,  we  see  that  it  is  impossible  for  hQ1  to  equal  0 
at  consecutive  times.  Thus,  there  are  only  three  possible  SI  states  which  we 
have  coded  as  follows  in  Figure  5: 

0  =  (0.  1)  .  1  =  (1.  0)  .  2  =  (1.  1) 

Since  hgj  is  a  deterministic  function  of  the  state  of  00  it  is  straightforward 
to  derive  the  way  in  which  Xg(n)  affects  the  transition  behavior  of  SI  (see 
Figure  5). 

As  in  all  of  our  models,  the  observation  z^(n)  is  associated  with 
t rans 1 1 1 ons  in  the  SI  subprocess  which  correspond  to  3-tuples  (hQj(n-l), 

h0i(n-3))  interactions  (corresponding  to  the  states  of  the  SI 
subprocess  before  and  after  the  transition).  Our  measurement  model  is  then 


S2  Model  (p=  0.00784) 

P 


transition  rates  are  as  given  in  Figure  1  with  x  =  13,  only 
if  the  S2  process  is  in  state  2.) 


the  set  of  conditional  probabilities 

Prfz^n)  |h01(n-l),h01(n-2),h01(n-3))  (4) 

Since  the  Viterbi  algorithm  provides  us  with  noncausal  estimates,  we  are  free 

to  build  some  noncausality  into  this  model.  Consequently,  we  have  chosen  to 

take  z^(n)  as  the  pass  1  estimate  at  time  n-2,  which  therefore  provides  an 

estimate  of  hQ^(n-2).  Thus  the  model  allows  us  to  capture  time  shifts  of  ±1.* 

The  specification  of  (4)  can  be  obtained  by  analysis  of  the  performance  of  the 

first  step  estimator.  We  have  estimated  these  quantities  by  simulating  that 

estimator  and  tabulating  the  results. 

We  now  must  consider  the  interactions  h^g(n)  initiated  by  Cl  and 

impinging  on  00,  i.e.  the  effect  of  V-wave  occurrences  on  CO.  There  is  a 

similarity  here  with  the  modeling  of  SO  in  the  first  pass  but  in  the  present 

context  we  also  have  the  estimates  from  pass  1  which  tell  us  something  about 

these  interactions.  Specifically,  since  we  used  the  exact  Cl  submodel  in  pass 

1.  we  can  deduce  estimates  of  h^  (see  (3)).  We  take  the^e  estimates  as  our 

observation  z^  for  pass  2  (without  any  augmentation  as  was  done  for  since 

the  first  step  estimator  used  an  exact  model  for  Cl  and  consequently  should 

produce  comparatively  accurate  estimates).  Also,  as  with  the  SI  submodel,  we 

need  to  model  possible  estimation  timing  errors,  so  again  we  take  the  state  of 

S2  to  be  a  set  of  the  most  recent  interactions,  in  this  case  (h^(n), 

2 

hio(n-l)).  In  this  example  it  is  impossible  for  h^Q  to  equal  1  at  two 

*Note  that  allowing  timing  errors  of  both  signs  is  needed  in  any  realistic 
algorithm  since,  for  example,  wave  detection  preprocessors  are  noncausal  in 
nature . 

^Note  that  there  is  some  asymmetry  in  comparison  with  the  SI  submodel  where 
the  state  was  lagged  one  step.  This  is  a  result  of  the  fact  that  in  the  SI 
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consecutive  times,  and  thus  we  can  code  the  feasible  S2  states  as 

0  =  (0,  0)  .  1  =  (0,  1)  .  2  =  (1,  0) 

In  this  example  the  CO  submodel  transition  probabilities  are  as  illustrated 
pictorial ly  in  Figure  1  for  Xg2(n)  =  0  or  1  and  incorporate  the  0.5 
probability  reset  to  state  25  when  Xg2(n)  =  2.  The  model  for  S2  is  also 
illustrated  in  Figure  5.  Note  that  as  with  SI,  there  is  a  parameter  p  to  be 
chosen  to  specify  the  S2  transition  probabilities.  This  parameter  was  also 
chosen  to  match  statistics  of  the  true  h^  process  (as  were  the  other 
probabilities  of  the  S2  process,  but  these  can  be  determined  trivially  from 
the  structure  of  the  Cl  submodel).  The  general  method  used  is  described  in 
the  next  section.  Finally,  the  observation  z2(n),  which  is  the  pass  1 
estimate  of  h^Q(n-l).  is  modeled  as  resulting  from  S2- transit ions.  Thus  again 
we  must  specify  a  distribution,  namely 

Pr(z2(n) |h10(n),  h1Q(n-l).  h1Q(n-2)) 
which  we  have  again  done  by  simulation. 

This  completes  the  specification  of  the  second  pass  model.  It  is  worth 
noting,  however,  that  there  are  several  things  that  have  been  left  out  of  the 
model.  The  most  glaring  of  these,  perhaps,  is  the  complete  absence  of  R-,  T-, 
and  V-waves.  For  the  pass  1  estimation  algorithm  we  argued  that  it  was 
reasonable  to  consider  omitting  P-waves  from  the  model  since  (a)  we  were 
focusing  most  attention  on  submodel  Cl  and  (b)  the  P-waves  were  of  low 
amplitude.  In  pass  2,  the  first  argument  holds  (here  we  are  focusing  on  CO), 

submodel,  h^^n)  is  a  deterministic  function  of  Xg(n).  Thus  for  the  state 

XQ(n)  to  correctly  "influence’’  the  next  transition  in  SI,  we  needed  to 

introduce  the  time  delay  in  defining  the  SI  state.  This  is  not  needed  in  S2, 
since  there  is  no  such  deterministic  coupling. 


but  the  latter  does  not,  since  R-,  T-,  and  V-waves  are  not  of  low  amplitude. 

In  the  general  procedure  described  in  the  next  section,  we  allow  for  the 
possibility  of  taking  such  waves  into  account  through  so-called  sub tractor 
submodels.  However,  as  the  results  in  this  section  and  in  [Doerschuk  1985] 
indicate,  for  ECG-type  models  such  as  the  one  considered  here  the  use  of  such 
models  is  of  doubtful  value.  Intuitively  the  reason  we  can  ignore  such  waves 
in  the  pass  2  estimation  algorithm  is  that  through  z^(n)  and  Zgfn)  we  are 
providing  indications  of  the  times  at  which  these  waves  occur.  Given  then  the 
coupling  between  these  waves  and  the  likely  times  of  P-waves,  captured  in  the 
original  00-C1  model  and  in  our  simplified  pass  2  00-S1-S2  model,  the  pass  2 
estimator  will  not  try  to  account  for  R-,  T-,  and  V-waves  by  placing  P-waves 
in  their  locations. 

A  second  issue  we  have  ignored  is  that  of  allowing  the  00  submodel  to 
influence  the  S2  submodel  (which  is.  in  our  approach,  autonomous).  Since  the 
CO  submodel  does  influence  the  Cl  submodel,  there  would  seem  to  be  some 
argument  for  doing  so.  However,  it  is  precisely  this  influence  that  is 
focused  upon  in  the  SI  submodel,  while  the  S2  submodel  focuses  on  that  part  of 
the  Cl  submodel,  dealing  with  V-waves,  which  is  unaffected  by  the  CO  submodel. 
Consequently,  while  our  general  modeling  methodology  allows  CO  to  influence 
S2,  it  is  not  necessary  to  include  this  bit  of  complexity  in  the  present 
context . 

Finally,  we  note  that  in  our  model  we  are  considering  and  Zg  to  be 
independent  measurements,  which  is  clearly  erroneous  since  they  are  both 
determined  by  the  pass  1  estimation  process.  One  can  certainly  construct  a 
somewhat  more  complex  model  involving  a  joint  distribution  of  z^,  z ^  given  the 
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combined  information  in  the  most  recent  transitions  of  SI  and  S2.  but  this  was 
not  found  to  be  necessary  here  (since  again  and  Zg  focus  on  different 
portions  of  the  overall  C0-C1  model). 

In  summary,  the  second  pass  of  our  estimation  procedure  consists  of  the 
minimum  probability-of-error  estimation  of  the  state  trajectory  of  the  model 
given  in  Figure  5  (using  the  Viterbi  algorithm)  given  the  EGG  measurement  and 
the  derived  measurements  z^  and  Zg  from  the  first  pass.  The  results  for  this 
example  are  given  in  the  third  row  of  annotations  in  Figure  2  where  we  have 
indicated  the  times  at  which  P-waves  were  estimated  to  have  occurred. 

Comparing  this  to  the  top  row  of  annotations  we  see  that  performance  is  quite 
good.  Note  that  the  erroneous  R.T-wave  pairs  from  pass  1  near  136.6  and  138.3 
seconds  did  not  lead  to  an  erroneous  P-waves  in  pass  2,  thanks  to  our  modeling 
of  Zj  which  incorporated  the  possibility  of  such  false  alarms.  Note  also  the 
occurrence  of  P-wave  timing  errors  (as  illustrated  near  80.2  and  99.9  seconds) 
all  of  which  underestimate  the  P-R  interval.  Finally,  note  that  it  is 
possible  in  our  model  (and  in  the  heart)  for  P-  and  V-waves  to  occur  nearly 
simultaneously  or  for  V-waves  to  preempt  an  already  occurring  P-wave  from 
initiating  a  normal  R-wave.  Having  knowledge  of  this,  the  pass  2  estimator 
will  attempt  to  insert  P-waves  when  the  timing  seems  likely  even  though  the 
presence  of  V-waves  may  obscure  the  P-wave  (indeed  this  is  precisely  how  a 
human  would  create  such  estimates).  An  example  of  correct  estimates  of  this 
type  can  be  found  near  99  seconds.  A  false  alarm  can  be  seen  near  82.6 
seconds,  and  a  missed  detection  near  83.3  seconds.  While  the  value  of  such 
estimates  is  suspect  (and  not  of  particular  consequence)  they  do  provide 
rather  graphic  examples  of  the  use  our  estimator  makes  of  the  timing  and 
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control  information  embedded  in  our  models. 

The  third  pass  of  the  estimation  process,  whose  purpose  is  to  provide 
Improved  and  consistent  estimates  of  ventricular  activity,  is  based  on  a 
model,  illustrated  in  Figure  6,  with  structure  analogous  to  that  of  pass  2 
(with  the  roles  of  submodels  CO  and  Cl  interchanged).  Specifically,  since  the 
focus  of  pass  3  is  on  submodel  Cl,  we  include  a  complete  version  of  this 
submodel.  Also,  from  pass  2  we  obtain  estimates  of  interactions  hjQ  initiated 
by  Cl  and  impinging  on  CO  (these  come  from  the  estimates  of  the  state  of  the 
S2  submodel  in  pass  2)  and  estimates  of  interactions  hQ1  initiated  by  00  and 
impinging  on  Cl  (these  come  directly  from  the  pass  2  estimate  of  the  state  of 
CO).  As  before,  since  only  an  aggregate  model  (S2)  is  used  in  pass  2  to 
estimate  h^Q,  we  augment  these  estimates  with  the  corresponding  estimate  of 
the  state  of  CO.  These  augmented  estimates  are  incorporated  in  pass  3  as 
observations,  denoted  by  Zj  in  Figure  6,  of  an  SI  submodel  whose  state  is  the 
most  recent  two  values  of  the  actual  interactions,  (h1Q(n-l).  h1Q(n-2)).  The 
pass  2  estimates  of  hg^  are  incorporated  in  pass  3  as  the  observations,  Zg,  of 
an  S2  submodel  whose  state  is  (h^^n),  h^^n-l)).  The  same  methods  as  in  pass 
2  were  used  to  compute  the  transition  probabilities  and  measurement 
distributions.  The  estimator  is  again  a  minimum  probability-of-error 
estimator  using  the  ECG  and  the  derived  measurements  z^,  Zg. 

The  result  of  applying  this  estimator  is  illustrated  in  the  fourth  row  of 
annotations  in  Figure  2.  The  final,  overall  estimate  consists  of  the  CO-state 
estimate  of  pass  2  (row  3)  and  the  Cl-state  estimate  of  pass  3  (row  4),  which 
are  combined  in  row  5.  Comparing  the  top  and  bottom  rows  we  see  that  the 
estimator  has  performed  quite  well.  Disregarding  the  initial  heartbeat  (which 
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was  missed  in  pass  3  because  of  the  specific  way  in  which  we  implemented  the 
initialization  of  the  latter  passes  of  our  algorithm)  all  R-,  T-,  and  V-waves 
were  detected  and  located  with  no  false  alarms.  Note  that  while  there  had 
been  several  false  R,  T-wave  estimates  in  pass  1.  these  have  been  completely 
eliminated  in  pass  3.  in  which  we  have  the  benefit  of  using  estimates  of 
OO-behavior  in  order  to  enforce  consistent  overall  estimation. 

The  estimation  of  P-wave  occurrences  is  also  quite  good.  Quantifying 
this  performance,  however,  is  an  interesting  question  itself,  since  one  is 
clearly  not  Just  Interested  in  estimation  errors  at  points  in  time  but  also  in 
timing  errors  at  {Mints  in  the  estimated  event  cycle  —  i.e.  an  estimation 
error  of  one  time  sample  in  locating  a  P-wave  should  not  be  thought  of  as  a 
missed  detection  but  rather  as  a  (not  particularly  troubling)  timing  error. 
Much  more  on  the  issue  of  performance  measures  for  event-oriented  estimation 
problems  can  be  found  in  [Doerschuk  1985].  This  example  does,  however, 
indicate  the  main  ideas.  In  examining  the  results  of  the  full  simulation  we 
find  that  there  are  only  two  Isolated  false  positive  P-wave  indications  and 
one  Isolated  false  negative  (neglecting  the  initial  heartbeat),  where  by 
"isolated"  we  mean  that  there  is  no  nearby  P-wave  in  the  true  or  estimated 
state  trajectories.  Given  that  there  are  230  heartbeats  in  this  simulation, 
these  correspond  to  a  false  positive  rate  of  .009  and  a  false  negative  rate  of 
.004.  There  are  also  23  other  paired  false  positives  and  negatives,  where  we 
have  used  the  criterion  of  associating  estimated  and  actual  P-wave  locations 
only  if  the  waveforms  at  these  locations  overlap.  This  corresponds  to  a 
paired  error  rate  of  0.10.  Note  that  in  our  model,  every  R-wave  must  be 
preceded  by  a  P-wave,  and  thus  this  pairing  is  to  be  expected  (we  discuss 


a  less  deterministically-coupled  model  in  [Doerschuk  1985]).  It  is  worth 
noting  that  in  each  of  these  paired  errors,  the  estimated  P-wave  location  was 
closer  to  the  R-wave  than  the  true  R-wave,  Indicating  a  bias  that  may  be 
removable  (and  is  most  likely  due  to  the  pass  2  estimator  correlating  the 
P-wave  with  the  initial  portion  of  the  R-wave). 

In  [Doerschuk  1985]  we  consider  a  variety  of  other  models  with  other 
types  of  variability  them  that  captured  in  the  model  of  Figure  1.  For 
example,  we  have  examined  models  with  transient  AV  block,  i.e.  models  in  which 
not  every  attempt  at  ventricular  excitation  leads  to  tin  R-wave  even  if  the 
ventricles  are  apparently  in  the  resting  state.  Because  of  the  additional 
freedom  in  the  model,  one  would  expect  some  drop  in  performance.  However  the 
drop  is  extremely  small  both  for  a  globally  optimal  estimator  and  for  the 
suboptimal  estimator  designed  based  on  the  principles  outlined  in  this  section 
and  formalized  in  the  next. 

3.  A  General  Design  Methodology 

The  example  described  in  detail  in  the  previous  section  illustrates  the 
major  elements  of  a  general  estimator  design  methodology  for  distributed 
Markov  chains.  In  this  section  we  describe  that  methodology.  Specifically, 
consider  the  estimation  of  an  interconnection  of  subprocesses,  which  we  denote 

CQ.  Cj . Cjj,  with  states  Xq.  Xj-.-.x^,  given  measurements  of  signals 

containing  signatures  corresponding  to  particular  state  transitions  in  these 
subprocesses.  The  interactions  among  these  subprocesses  are  defined  as 
follows.  Let  hAj(n)  denote  the  interaction  initiated  by  Ci  and  impinging  on 
Cj  at  time  n.  This  interaction  is  a  deterministic  function  of  x^(n),  and  the 


transition  probabilities  of  Cj  are  deterministic  functions  of  (h^fn)  |i^j} . 

The  assumption  is  that  the  set  of  possible  transition  probabilities  for  each 
Cj  (and  thus  the  set  of  possible  values  of  (h^ j(n) |i/j})  is  quite  small. 

Our  overall  estimator  consists  of  an  interconnection  of  local  estimators 
(LE's).  each  of  which  focuses  on  the  estimation  of  one  of  the  subprocesses. 
Because  of  the  existence  of  interactions  with  and  events  in  the  observed  data 
due  to  other  subprocesses,  each  LE  not  only  must  take  these  effects  into 
account  in  its  model  but  also  must  communicate  with  the  other  LE's  in  order  to 
determine  a  consistent  and  accurate  set  of  estimates.  Consider  first  the 
initial  pass  through  the  data,  at  which  time  the  LE's  have  no  previous 
information  to  communicate.  Let  us  focus  on  the  LE  for  a  specific  submodel, 
namely  Cj .  In  order  for  this  LE  to  construct  its  initial  estimate  it  in 
general  it  will  need: 

(a)  A  complete  model  of  the  subprocess  Cj  on  which  it  is  focused. 

(b)  A  model  of  the  interactions  impinging  on  C  . 

J 

(c)  A  model  of  the  waveforms  generated  by  the  other  submodels. 

The  model  referred  to  in  (b)  is  called  an  SO  submodel,  and  a  major 

objective  is  to  make  it  as  simple  as  possible  in  order  to  keep  the  LE  as 

3 

simple  as  possible.  In  the  example  in  the  previous  section  and  throughout 

^There  are  two  distinct  ways  in  which  one  can  perform  this  modeling  step  and 
several  that  follow.  In  particular,  in  this  section  we  describe  the 
construction  of  a  single  SO  submodel  capturing  the  interactions  impinging  on 
Cj  from  all  other  subporcesses.  In  [Doerschuk  1985]  an  analogous  approach  is 

described  for  constructing  separate  SO  submodels  for  the  interactions 
initiated  by  each  of  the  other  subprocesses. 


our  work  we  have  taken  the  states  of  the  SO  submodel  to  be  In  one-to-one 
correspondence  with  the  possible  values  of  the  N-tuple  {b^(n)|i=J).  Whether 
this  or  a  more  complex  model  is  used  to  describe  the  dynamics  of  these 
interactions,  the  question  remains  of  choosing  the  transition  probabilities 
for  the  90  submodel.  Our  approach  has  for  the  most  part  been  to  match  these 
one-step  translation  probabilities  to  the  actual  steady-state  versions  within 
the  original  process,  that  Is.  to 

lim  Pr[{h1J(n)|i#j}|{hiJ(n-l)|i*j}.  (hj^n-l)  =  hJt  |i/j}]  (5) 

This  computation  deserves  some  comment.  Note  first  that  we  have  Included 
conditioning  on  (hj^n-l)  |i*J} .  which  reflects  the  Influence  has  on  the 
other  subprocesses.  This  results  in  the  transition  probabilities  of  SO  being 
Influenced  by  the  state  of  Cj.  Again  we  typically  expect  this  Influence  to 
manifest  Itself  as  a  small  number  of  possible  values  for  a  small  subset  of  the 
transition  probabllles  (e.g.  in  our  case  study  only  the  parameter  p  in 
Figure  4  is  influenced,  and  it  only  takes  on  two  values).  Note  that 
conditioned  on  {hji(n).|i#j} .  {x^njji^J}  is  a  Markov  chain.  However  this  is 
not  typically  the  case  for  the  highly  aggregated  (h^n)  | i*j} .  so  that  the 
limit  in  (5)  is  not  a  completely  trivial  computation.  On  the  other  hand,  once 
we  compute  the  ergodic  probabilities  for  {x^njji^j},  it  is  straightforward  to 
compute  (5) . 

Typically  for  models  with  infrequent  changes  in  interactions,  most  of  the 
transition  probabilities  specified  in  (5)  are  0  or  1,  and  there  are  only  a  few 

A 

parameters  (such  as  p  in  Figure  4)  for  which  this  computation  is  necessary. 

4 

Indeed  for  al 1  of  the  cases  considered  in  [Doerschuk  1985],  the  model  was 
exactly  as  in  Figure  4  (with  different  values  of  p),  since  in  all  of  our  cases 
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Finally,  we  note  that  there  are  some  cases  in  which  the  matching  of  the 
steady-state  statistic  (5)  may  be  inappropriate.  In  particular,  in  using  (5) 
we  are  in  essence  assuming  that  the  transition  probabilities  of  {x^(n)|i*j}  do 
not  change  very  frequently  (so  that  steady-state  is  actually  achieved)  —  i.e. 
that  the  time  variations  observed  in  the  actual  xj(n)  process  do  not  lead  to 
frequent  changes  in  the  interactions  h^(n)  (which  are  assumed  to  be 
constant  in  (5)).  We  refer  the  reader  to  [Doerschuk  1985]  for  examples 
violating  this  assumption  and  in  which  we  must  set  the  SO  transition 
probabilities  in  a  different  manner.  Note  that  this  assumption  is  in  fact 
violated  in  our  case  study.  In  particular,  while  it  is  certainly  true  that 
hjQ  =  0  for  long  periods  of  time,  since  this  corresponds  to  Xj  being  any  state 
other  than  13.  h^  =  1,  corresponding  to  x^  =  13.  cannot  possibly  occur  at  any 
two  consecutive  times.  In  this  case,  since  h^  =  1  corresponds  to  a  reset  of 
CO  to  state  25.  and  since  all  states  in  00  other  than  0  correspond  to  hQ1  =  1, 
it  is  reasonable  to  reset  the  state  of  SO  to  1  whenever  x^^  =  13.  This  is  what 
is  specified  in  Figure  1  and  what  we  would  calculate  from  (5).  Thus  (5)  is 
often  useful  even  if  the  assumption  on  which  it  is  based  is  violated. 

Let  us  now  discuss  the  model  referred  to  in  (c).  This  is  one  of  the 
subtractor  submode 1 s .  denoted  by  S3,  referred  to  in  the  previous  section, 
which  is  incorporated  in  order  to  keep  the  LE  from  interpreting  waveforms 


generated  by  other  submodels  as  coming  from  C 


In  essence  what  we  would  like 


to  do  is  to  present  the  LE  with  observations  containing  only  those  signatures 


generated  by 


Since  this  is  not  possible,  we  equip  the  LE  with  a  mechanism 


there  have  been  only  two  interaction  values,  one  of  which  could  not  occur  at 
consecutive  times. 
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for  estimating  when  other  signatures  have  occured  so  that  it  can  in  effect 
subtract  out  their  effects.  In  general,  one  can  construct  a  separate  S3 
subendel  for  each  signature  not  initiated  by  Cj.  While  it  would  certainly  be 
possible  to  couple  these  subprocesses  with  the  Cj  and  SO  submodels,  we  have 
obtained  good  results  with  a  simpler  structure  in  which  each  S3  submodel  is  a 
completely  autonomous,  aggregated  process  that  produces  interarrival 
statistics  for  the  wave  of  interest  identical  to  those  produced  by  the  exact 
model.  More  precisely,  we  have  always  taken  S3  submodel  to  be  of  the  form 
shown  in  Figure  7.  The  two  parameters  are  chosen  as  follows.  Let  Tgg(n) 
denote  the  time  between  the  nth  and  ( n+ 1 ) s t  occurrence  of  the  signature  S  in 
the  original  process.  Then  we  choose  p  and  q  to  match  the  probability  that 
signatures  occur  at  successive  times  and  the  mean  time  between  successive 
signatures.  That  is 


(6) 

(7) 


probabilities  of  the  full  model.  In  most  cases  Pr[TOQ(n)  =  1]  =  0,  so  that 


lim  EOggfn)]  -  1 
n-*» 


Therefore,  in  our  general  design  methodology  we  construct  each  initial  LE 


model  using  C  .  SO.  and  S3  components  as  illustrated  in  Figure  8  and  compute 

J 


the  initial  pass  minimum  probabi 1 i ty-of-error  estimates  for  each  LE.  We  are 


then  in  a  position  to  consider  a  ref inement  pass,  in  which  each  LE  reprocesses 


the  data,  together  with  information  provided  from  the  initial  passes  of  the 


LE's. 


waveforms  initiated  by  Cj 


waveforms 
initiated  by 
other  submodels 


The  structure  of  a  general  LE  model  for  an  initial  pass 


The  LE  for  subprocess  Cj  will  in  general  need  the  following  elements  in 
its  model  for  a  refinement  pass: 


(1) 


A  complete  model  of  C 


j' 


(2) 


A  model  of  the  information  provided  by  the  previous  pass 
concerning  interactions  initiated  by  Cj. 


(3)  A  model  of  the  information  provided  by  the  previous  pass 
concerning  interactions  impinging  on  C^. 


(4)  A  model  of  the  information  provided  by  the  previous  pass 

concerning  times  of  occurrence  of  waveforms  generated  by  the 
other  subprocesses. 


Note  that  (2)  and  (3)  together  corresponds  to  (b)  in  the  initial  pass.  We 

have  chosen  to  split  these  two  models  here  as  (i)  it  makes  the  modeling  of  the 

information  simpler  and  (ii)  the  information  referred  to  in  (2)  and  (3)  may 
typically  come  from  different  sources  or  be  of  very  different  accuracy  or 
structure  (given  that  each  LE  has  an  accurate  model  of  its  own  subprocess  but 
only  highly  aggregated  models  of  the  others). 

As  discussed  in  the  previous  section,  the  models  referred  to  in  (2)  and 
(3),  which  we  refer  to  as  SI  and  S2  submodels,  respectively,  must  capture  the 
timing  and  estimation  uncertainties  from  the  previous  pass.  Each  accomplishes 
this  by  taking  as  its  state  space  a  moving  window  of  the  most  recent 

interactions.  In  particular,  the  state  of  the  SI  submodel  consists  of  a 

window  of  the  most  recent  values  of  the  N-tuple 

{hjji^j} 

while  the  state  of  the  S2  submodel  is  a  window  of  the  most  recent  values  of 
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the  N- tuple3 
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An  objective  in  designing  these  models  is  to'  keep  the  window  lengths  and 
small  in  order  to  minimize  state  space  size.  This  desire  is  balanced  by  the 
need  to  model  estimation  timing  errors  (since  the  maximum  such  symmetric  error 
that  can  be  modeled  corresponds  to  1/2  the  window  length).  In  our  work  we 
have  always  taken  this  window  length  equal  to  2. 

Consider  next  the  dynamics  of  the  SI  and  S2  submodels,  both  of  which  are 
specified  in  a  straightforward  manner.  In  particular,  since  each  h^n)  is  a 
deterministic  function  of  xj(n)  and  since  the  full  Cj  model  is  used  by  the  LE. 
the  SI  dynamics  are  essentially  a  shift  register  memory.  That  is,  given  Xj(n), 
the  transition 

Xgjfn)  =  (hjj(m)  ( i*j ,  n^n-I^ . n-1} 

1 

Xg^n+l)  =  (hji(m)  | i*j  ,  m=n-K2+l, - n) 

is  deterministic  (i.e.  for  each  present  state  there  is  one  next  state  [whose 
identity  depends  on  x^(n)]  that  SI  will  occupy  with  probability  1). 

The  dynamics  of  the  S2  submodel  are  not  quite  this  trivial. 

Specifically,  as  in  the  SO  submodel,  we  choose  the  transition  probabilities  to 
match  those  in  the  original  process.  In  particular,  in  analogy  with  (5),  we 
can  choose  these  to  equal 

lim  Prtlhjjfm)  |i*j,  m=n-K2+2 . n+1}  |{hij(m)  |  i*j  ,m=n-K2+l - ^.{h^fn)  |i*j}] 


Recall  from  the  previous  section  that  there  is  some  asymmetry  in  the  windows 
here,  with  the  window  for  SI  stopping  at  time  n-1,  and  the  window  for  S2 
stopping  at  time  n. 


m 
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By  including  the  conditioning  on  {h^fnjli/j}  we  can  capture  the  interactions 

initiated  by  Cj  and  impinging  on  the  other  subprocesses  (and  therefore,  in  the 

LE  model,  on  S2).  However,  as  discussed  in  the  previous  section,  the  effects 

of  these  interactions  are  the  primary  concern  of  the  SI  submodel,  and  thus  it 

is  worth  seeking  and  typically  possible  to  find  a  far  simpler  model.  In  fact 

throughout  our  work  we  have  been  able  to  completely  eliminate  the  influence  of 

Cj  on  S2  (which  then  operates  autonomously,  generating  the  interactions  that 

impinge  on  C^).  This  can  be  done  by  using  (9)  with  {hjj(n)|i^j}  set  equal  to 

6 

the  values  that  represent  the  most  usual  interaction.  Another  approach  is  to 
compute  the  average  of  (9)  over  the  possible  values  of  {hjj(n)|i*j}  using 
their  ergodic  probabilities.  In  our  work  we  have  used  the  latter  of  these  two 
methods. 

Consider  next  the  modeling  of  the  ’’measurements”  provided  by  the  previous 
data  pass.  With  respect  to  SI,  we  have  in  general  the  following  sources  of 
information  concerning  the  interactions  initiated  by  C^: 


(i)  The  previous  state  estimate  of  from  its  associated  LE.  From 
this  we  can  directly  compute  an  estimate  of  (h^fn)  |i?fj} . 


(ii)  The  augmented  interactions  from  each  of  the  other  LE’s.  These 
consist  of  the  estimate  of  the  interaction  impinging  on  the  CL 

submodel  associated  with  each  LE  (obtained  from  the  aggregated  SO 
submodel  used  by  the  LE)  and  the  corresponding  C^-state  estimate 

(see  the  previous  section  for  a  discussion  of  why  this  augmenting 
information  is  included). 


In  our  ECG  examples  this  corresponds  to  no  attempt  at  interprocess 
excitation,  as  such  electrical  excitations  occur  over  relatively  short  time 
periods  (usually  a  single  time  sample). 
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Together  this  information  forms  a  measurement,  which  we  denote  z^(n),  and, 
since  we  associate  measurements  with  transitions,  we  model  the  information 
contained  in  Zj(n)  by 

Pr[Zj(n) |{h  (m)i*J,  msn-Kg-l . n-1}]  (10) 

As  discussed  in  the  previous  section,  since  we  have  introduced  some 
noncausality  with  the  Viterbi  algorithm,  we  have  the  flexibility  of  doing  so 
here  (and  the  need  in  order  to  model  positive  and  negative  timing  errors). 
Consequently,  we  take  z^(n)  to  be  the  previous  pass  estimates  indicated  in  (i) 
and  (ii)  evaluated  at  time  n-l-^^.  Finally,  while  it  is  possible  to  devise 
analytical  methods  to  obtain  approximations  for  (10).  we  have  found  it  easier 
to  evaluate  these  distributions  by  simulation. 

For  S2,  we  have  the  following  sources  of  information  concerning 
interactions  impinging  on  Cj= 


(i)  The  augmented  estimated  interaction  provided  by  the  previous  pass 
of  the  LE  for  C^. 


(ii)  The  estimated  state  of  each  provided  by  the  associated  LE. 
From  these  we  can  directly  compute  estimates  of  each  h^n). 


This  information  forms  the  measurement  z2(n),  which  is  modeled  via 

Pr[z2(n) | (hA j(n) j i^j ,  i=n-K2 - ,n}]  (11) 

As  in  the  case  of  SI,  we  introduce  some  noncausality  by  taking  z2(n)  to  be  the 
previous  pass  information  evaluated  at  n  -  ^^2,  and  we  determine  (11)  by 
simulation. 

Finally,  consider  modeling  the  information  available  from  the  previous 
pass  concerning  waveforms  generated  by  other  submodels.  Each  such  waveform  is 
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modeled  by  what  we  call  an  S4  submodel,  which  is  a  second  class  of  subtractor 


models  used  in  refinement  passes.  S4  submodels  resemble  S2  submodels  in 


structure  and  principle.  Consider  an  S4  submodel  corresponding  to  a 


particular  waveform  generated  by  submodel  C^.  The  measurement  z^(n)  provided 


by  the  previous  pass  LE  for  is  a  sequence  of  binary  annotations  —  0  if  the 


LE  estimates  that  the  particular  waveform  was  not  initiated  at  that  time 


sample  and  1  if  the  estimate  is  that  the  waveform  was  generated.  The  state  of 


the  S4  submodel  is  a  window  of  the  most  recent  true  values  of  these  binary 


annotations.  As  with  S2.  the  transition  rates  of  this  model  are  chosen  to 


match  the  corresponding  transition  rates  of  sequences  of  binary  annotations  in 


the  full  model.  If  the  counterpart  to  (9)  is  used,  the  S4  model  will  in 


general  be  influenced  by  C. .  Again  as  in  the  case  of  S2,  we  have  typically 


simplified  this  model  so  that  S4  is  autonomous,  by  averaging  out  the 


Cj-dependence  using  the  ergodic  distribution  for  Xj(n). 


The  output  of  the  S4  chain  is  a  sequence  of  occurrences  of  the  waveform 


being  modeled.  Such  outputs  occur  at  all  S4  transitions  to  states  with  a  1  as 


the  most  recent  annotation.  The  auxiliary  observation  z^(n)  is  again  modeled 


via  a  probability  distribution  conditioned  on  the  most  recent  S4  transition. 


We  have  determined  distributions  of  this  type  via  simulation. 


The  structure  of  the  models  on  which  each  LE  refinement  pass  is  based  is 


depicted  in  Figure  9.  Again  the  Viterbi  algorithm  is  used.  In  general  one 


can  envision  making  several  refinement  passes,  with  the  final  estimate 


consisting  of  the  collection  of  C^-state  estimates  from  the  final  passes  of 


the  corresponding  LE's.  The  primary  purpose  of  the  refinement  passes  is  to 


improve  the  accuracy  and  consistency  of  this  set  of  estimates.  In  particular, 
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if  one  implemented  a  single,  optimal  estimator  for  the  full  process,  one  would 
know  for  certain  that  all  transitions  present  in  the  final  state  estimate 
would  be  consistent  (i.e.  have  nonzero  probability  in  the  full  process).  When 
one  uses  a  collection  of  distributed,  simpler  LE’s.  there  is  no  such 
guarantee,  but  the  coordination  made  possbile  by  refinement  passes  makes  the 
occurrence  of  inconsistent  estimates  far  less  likely.  In  fact,  in  our  studies 
consistency  has  not  proven  to  be  a  problem. 

The  complete  procedure  we  have  described  requires  the  implementation  of  a 
full  set  of  LE*s  for  the  initial  pass  (based  on  models  as  in  Figure  8)  and 
subsequent  refinement  passes  (each  based  on  its  own  model  as  in  Figure  9).  As 
in  our  example  in  Section  2,  it  is  typically  possible  to  simplify  this  design 
considerably.  First  of  all.  as  in  our  example,  for  each  LE  it  is  often  not 
necessary  in  the  initial  pass  to  include  subtractor  submodels  S3  for  waveforms 
of  low  SNR  compared  to  the  waveforms  generated  by  the  submodel  corresponding 
to  the  LE.  Also,  as  we  showed,  it  may  not  be  necessary  to  include  any  S4 
submodels,  since  the  information  provided  through  SI  and  S2  submodels 
essentially  provides  timing  information  that  allows  the  LE  to  avoid  intervals 
in  which  these  interfering  signatures  may  appear.  Doerschuk  [1985]  presents 
comparative  results  with  and  without  S3  and  S4  submodels  that  support  these 
simplif ications. 

It  is  also  typically  possible  to  eliminate  many  of  the  LE’s  from  each 
pass.  For  example,  in  the  initial  pass,  one  typically  would  implement  LE's 
only  for  submodels  generating  the  higher  SNR  signatures  (such  as  R-waves),  as 
the  performance  of  initial  pass  LE’s  for  other  submodels  with  only  low  SNR 
signatures  (or  no  signatures,  as  is  the  case  for  some  rhythm  models  described 
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in  [Doerschuk  1985,  1986])  will  generally  be  unsatisfactory.  Also,  in  order 
to  achieve  consistency,  we  do  not  need  to  refine  all  LE’s  in  subsequent 
passes.  In  particular  we  typically  can  implement  an  alternating  iterative 
structure  much  as  in  the  example  in  which  we  initially  estimate  the  C  with 

J 

high  SNR  signatures,  then  use  these  estimates  to  assist  in  estimating  only  the 
remaining  during  the  next  pass;  these  estimates  can  then  be  used  in  turn 
during  the  following  pass  in  the  re-estimation  of  the  from  the  initial  pass 
in  order  to  improve  the  accuracy  and  consistency  of  the  C^-state  estimates. 
Note  that  in  addition  to  eliminating  entire  passes  of  LE’s,  such  a  structure 
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reduces  the  quantity  of  z^  and  z ^  measurements  to  be  processed  by  the 
remaining  LE’s.  In  fact,  the  full  set  of  such  information  described 


previously  has  some  redundancy,  reflecting  the  fact  that  perhaps  not  all  of 


this  intermediate  processing  is  needed.  The  structure  described  above 
simplifies  the  design  by  removing  these  redundant  sources  of  information. 
Doerschuk  (1985)  presents  results  favorably  comparing  reduced  designs  of  this 
type  to  estimators  incorporating  more  or  all  of  the  LE’s  at  each  stage. 


4.  Conclusions 

In  this  paper  we  have  presented  a  methodology  for  the  distributed 
estimation  of  interconnected  finite-state  processes  given  the  observation  of 
signals  containing  waveforms  initiated  by  events  in  the  various  processes. 

Our  motivation  for  examining  this  class  of  models  is  the  problem  of  automated 
EGG  analysis,  but  the  methods  and  concepts  we  have  developed  are  of  potential 
use  in  a  variety  of  other  applications  such  as  the  monitoring  of  distributed 
power  networks. 


no 


The  approach  we  have  developed  highlights  the  major  issues  that  must  be 
addressed  in  designing  distributed  estimators,  namely  the  aggregated  modeling 
of  the  interactions  between  other  portions  of  the  overall  process  and  the 
particular  subprocess  being  estimated  and  the  dynamic  modeling  of  the 
information  provided  by  other  estimators  as  part  of  the  process  of  producing 
coordinated,  consistent  estimates  of  all  the  subprocesses.  We  have  presented 
systematic  procedures  for  constructing  these  models  that  can  in  fact  be  used 
as  the  basis  for  a  completely  automated  estimator  design  procedure  [Doerschuk 
1985]. 

In  order  to  illustrate  the  various  elements  of  our  design  process  we  have 
presented  a  case  study  corresponding  to  the  tracking  of  a  particular  cardiac 
rhythm  using  synthetic  data.  The  results  presented  indicate  the  potential  of 
this  design  method.  Two  major  issues  remain  to  be  considered,  however,  before 
a  complete  EOG  rhythm  analysis  system  can  be  constructed.  In  particular, 
while  our  distributed  design  yields  estimators  with  far  more  modest 
computational  demands  than  the  corresponding  optimal  estimator,  several  steps 
can  be  taken  to  simplify  these  computations  even  more.  First,  as  mentioned 
previously,  it  is  possible  to  construct  pipelined  versions  of  our  multi-pass 
estimators  in  which  all  passes  are  performed  at  the  same  time  rather  than  in 
sequence.  This  achieves  a  several-fold  increase  in  processing  throughput. 
Also,  the  nature  of  the  models  arising  in  EOG  analysis  offer  another 
possibility  for  simplification.  Specifically,  these  finite-state  processes 
typically  display  multiple  time  scale  behavior  (as  actual  signature-initiating 
events  occur  at  a  far  lower  rate  than  the  sampling  rate  needed  to  capture 
interprocess  timing).  Consequently,  it  may  be  possible  to  use  results  on 
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hierarchical  aggregation  of  processes  with  several  time  scales  [Coderch  1983] 
to  construct  more  efficient  estimators  that  not  only  display  the  spatial 
decomposition  we  have  exploited  here  but  also  the  time  scale  decomposition  of 
these  processes. 

Finally,  it  is  important  to  realize  that  the  problem  of  rhythm  tracking 
addressed  here  is  only  a  first  step  in  a  rhythm  diagnosis  sytems. 

Specifically  in  such  a  system  one  wishes  to  identify  the  underlying 
distributed  process  model  from  a  set  of  such  models  representing  different 
cardiac  rhythms.  As  in  standard  system  identification  problems,  the 
computation  of  the  likelihoods  for  a  set  of  models  can  be  performed 
efficiently  using  the  estimates  produced  by  estimators  based  on  each  of  the 
models  (e.g.  see  [Gustafson  1978a,  b]  for  an  application  of  this  idea  to  EOG 
rhythm  analysis  based  on  R-wave  location  data  only).  Doerschuk  (1985) 
describes  an  approach  to  constructing  such  likelihoods  based  on  the  outputs  of 
a  set  of  estimators  of  the  type  described  in  this  paper,  but  work  remains  to 
be  performed  to  test  this  method  and  to  develop  efficient  implementations. 
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